function [ vh ] = vcycle( vh, fh, v1, v2, current, lastGrid )
%VCYCLE Summary of this function goes here
%   Detailed explanation goes here

%global h;
global A;
global P;
%global dimm; % for debug

%vh
%fh
%A{current}

%A = prepareA(floor(sqrt(length(fh))), floor(sqrt(length(fh))), h);

% Step 1: Relax v1 times on Au=f
for i = 1 : v1
    vh=jacobi(vh, fh, A{current});
end

% Step 2: If current iteration is not the last one:
if current < lastGrid
    %f2h = rest(fh-A{current}*vh);
    %fh-A{current}*vh
    %P{current}'*(fh-A{current}*vh)
    %(0.5 .* P{current}')*(fh-A{current}*vh)
    f2h = P{current}'*(fh-A{current}*vh);
    v2h = zeros(length(f2h), 1);
    v2h = vcycle(v2h, f2h, v1, v2, current+1, lastGrid);
    % Step 3: Correct vh:
    %vh = vh + prolong(v2h);  
    %v2h
    %P{current}*v2h
    vh = vh + P{current}*v2h;   
end

% Step 4: Relax v2 times on Au=f
for i = 1 : v2
    vh=jacobi(vh, fh, A{current});
end
%figure; mesh(reshape(vh,dimm{current}));
end
